function eq = solve_eq_ss_bisect(param,glob,options,norm)

% Finds initial steady state
% With consumption normalized to 1

D1 = .001;
D2 = 5;

% Impose that C = 1. 
param.C = 1;   

Ds    = nan(options.itermax_DC,1);
% 
% if nargin == 5
%     param.D = eq0.D_out;
%     options.cresult = eq0.c;
% else
%     param.D = 19;
% end
%     
%     
 for t = 1:options.itermax_DC
     
    param.D = (D1 + D2)/2;
    
    Ds(t) = param.D;        
    
    param.w = param.omega*param.C; % HH FOC

    if norm == 1
        param.mu_f = .9*param.mu_entry;
    end
    
    tic
    eq = solve_cL(param,glob,options);     % solve incumbent & entrant pbms        
    eq = setup_transition_matrices(eq, param,glob,options); % set up transition matrices
    
    if norm == 1
        c2   = reshape(eq.v.v2, glob.n(1), glob.n(2));
        v2   = griddedInterpolant(glob.B, glob.S,c2, 'linear');   
        ve   = valfunc_E(v2, glob.agrid,param,glob,options); 
        tosolve = @(mu) .5 - logncdf(ve(floor(length(glob.agrid)/2)), mu, param.s_entry);        
        param.mu_entry = bisect(tosolve, -100, 100);
    end
    
    disp(param.mu_entry)

    
    eq = find_stationary_dist(param.M, eq, param,glob,options); % find stationary dist
    eq = find_agg(eq, param, glob, options); 
    toc;

    
    options.cresult = eq.c;

    subplot(2,1,1)
    hold off
    plot(Ds); title('Guesses of D'); if (t > 5); end;
    subplot(2,1,2)
    plot(glob.bgridf, eq.Lp); title('Distribution of customer base'); drawnow
    drawnow;
    
   diff = max(abs(param.D - eq.D_out)/param.D);
   if diff < 1e-6
        break
   end

   Douts(t) = eq.D_out;
   
   subplot(2,1,1)
   hold off
    line([0, t], [D1, D1], 'Color', 'Red')
    line([0, t], [D2, D2], 'Color', 'Red')
   hold on;
   plot(1:t, Douts, 'kx', 'MarkerSize', 5); drawnow;
    

  if eq.D_out < param.D
      D2 = param.D; 
      %D1 = max(eq.D_out, D1);
   else
      D1 = param.D;
      %D2 = min(eq.D_out, D2);
  end
 
   format short
   fprintf('Diff was %0.3e, new D is %.3f \n', diff, param.D)
   fprintf('D_out was %.3f \n', eq.D_out)
   fprintf('Exit rate was %.3f \n', eq.p_exit'*eq.L/sum(eq.L))
   fprintf('Entry rate was %.3f \n', sum(eq.G)/sum(eq.L))
   
end
    disp('hi')

    eq.mu_f     = param.mu_f;
    eq.mu_entry = param.mu_entry;

    eq.param = param;
    eq.glob = glob;
    eq.options = options;
    
end